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ABSTRACT 

The  nonlinear  deflections  of  a  thin  elastic  simply- 
supported  rectangular  plate  are  studied.   The  plate  is 
deformed  by  a  compressive  thrust  applied  along  the  short 
edges.   For  the  boundary  value  problem  considered  we  prove 
that  the  plate  cannot  buckle  for  thrusts  less  than  or  equal 
to  the  lowest  eigenvalue  of  the  linearized  buckling  problem. 
For  larger  thrusts  approximate  solutions  of  the  von  Karman 
equations  are  obtained  by  an  accelerated  iteration  method. 
Each  iterate  is  numerically  evaluated  by  a  finite  difference 
procedure.   Using  this  method  approximate  solutions  are  ob- 
tained for  thrusts  considerably  larger  than  the  lowest 
eigenvalue.   These  solutions  bifurcate  from  the  eigenvalues 
of  the  linearized  problem.   In  addition,  an  asymmetric 
solution  is  found  which  appears  to  branch  from  a  previously 
bifurcated  solution. 

The  extensive  numerical  results  are  used  to  study  the 
formation  of  boundary  layers  and  the  related  problem  of  the 
plate's  ultimate  load.   On  the  basis  of  the  numerical  results, 
an  energy  mechanism  is  proposed  to  explain  a  "mode -Jumping" 
phenomenon  which  has  been  previously  observed  in  experiments. 
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In  Sec.  3  of  this  paper,  a  part  of  the  above  conjecture 
is  proved.   We  show  that  there  are  no  buc!>.led  states  for 
thrusts  less  than  or  equal  to  the  lovrest  eigenvalue.   In 
this  analysis  the  nonlinear  von  Karman  plate  equations  [4] 
are  employed.   On  the  boundary  of  the  simply  supported  plate, 
the  tangential  "membrane"  displacements  are  taken  as  pro- 
portional to  the  tangential  coordinate,  and  the  normal 
"membrane"  stresses  on  the  unloaded  edges  are  assumed  to 
vanish. 

For  larger  thrusts,  approximate  solutions  of  the  von 
Karman  equations  are  obtained  by  an  accelerated  iteration 
procedure  which  is  related  to  one  previously  employed  in 
studies  of  the  nonlinear  axisymmetric  bending  and  buckling 
of  circular  plates  [5]  and  spherical  caps  [6],   Numerical 
approximations  to  the  bifurcating  solutions  are  obtained. 
In  addition,  vje  have  determined  another  buckled  state,  which 
appears  to  branch  from  a  bifurcating  solution,  and  which  we 
call  the  asymmetric  mode.   These  non-unique  solutions,  for  a 
wide  range  of  thrusts  considerably  in  excess  of  the  buckling 
load,  permit  a  study  of  the  development  of  boundary  layers. 

It  is  of  interest  to  determine,  for  a  fixed  thrust, 
which  of  the  many  buckled  states  a  "real"  plate  prefers. 
Experimental  results  (see  e.g.  [?])  show  that  there  exist  edge 
t^hrusts  at  which  plates  suddenly  "jump''  from  one  state  to  an- 
other. 
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This  indicates  that  the  preferred  state  m?,;/  vary  with  the 
thr^ist.   On  the  bo.sis  of  the  num-frrical  results,  v;e  propose 
in  Z-2C.    6  an  energy  mcC'::inlfent  to  describe  the  ^v:''^jlnn:  . 
phenomenon  and  we  sug£;ec;t  lower  Vjound's  on  the  ed.^^e  thrust 
at  vJhlch  it  can  occur.   The  pvechcnism  is  the  san-i  as  the  one 
pr'r:vic\:,ily  proposed  to  explain  the  snop  buclcling  of  spherical 
caps  [8] . 

With  suitable  modifications  the  present  numerical 
method  can  be  applied  to  other  flat  plate  buckling  problems. 
With  a  slight  modification  of  our  procedure  numerical 
solutions  for  the  nonlinear  bending  of  rectangular  plates 
due  to  the  uniform  lateral  pressure  and  the  nonlinear  buckling 
of  axially  compressed  cylindrical  shells  have  been  obtained. 
These  results  will  be  reported  elsewhere. 

In  previous  investigations  of  rectangular  plate  buckling 
other  boundary  value  problems  for  the  von  Karman  equation 
have  been  considered.   Perturbation  methods,  energy  methods 
and  series  expansions  were  used  to  obtain  approximate  solutions. 
For  references  and  discussions  see  the  texts  by  Tlmoshenko  [1] 
and  Bleich  [9]  and  the  papers  [2],  [7],  [10]  and  [11]. 
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2.   Formulation. 

Let  the  middle  surface  of  the  plate  coincide  with  the 
X,  Y  plane  so  that  the  undeformed  plate  occupies  the  region: 
|Z|  ^  t/2  ,  O^X^a,  0<Ylb.   The  plate  is  deformed  by 
a  constant  compressive  thrust  of  magnitude  T,  applied  normal 
to  the  edges  X  =  0,  a.   The  faces,  Z  =  +  t/2,  are  talten  as 
force  free . 

The  deformations  of  the  plate  are  assumed  to  be  adequately 
described  by  the  nonlinear  von  Karman  equations  [4].   These 
are  two  coupled  fourth  order  partial  differential  equations 
for  the  dependent  variables  W(X,Y)  and  P(X,Y).   W(X,Y)  is 
the  Z  displacement  of  the  plate's  mldplane  and  F(X,Y)  is  a 
stress  function.  The  von  Karman  equations  are  given  by, 

_2 

tA  W  =  P,yyW>xX  ^  ^'XX^^'YY  "  ^^'XY^^'XY  ' 


(2.1) 


_2       2 

A  F  =  E[M,^-W,j^^W,yy]  > 


where  E  and  v  are  Young's  mo  lulus  and  Poisson's  ratio, 
A  is  the  Laplacian  with  respect  to  X  and  2',  and 
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A  comma  denotes  partial  differentiation  with  respect  to  the 
subscripted  variables  that  follow.   The  "membrane"  stresses 
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cT  (X,Y),  6   (X,Y)  and  Q'     (X,Y)  are  expressed  In  terms  of 
X        y  xy 

F(X,Y)  by  the  equations, 

(2.2)     (S^   =  P,YY^  "^y  =  P>xX'   "xy  ==  "^'XY    ' 

Let  U(X,Y)  and  V(X,Y)  denote  the  midplane  displacements  in  the 
X  and  Y  directions.  They  are  related  to  the  membrane  stresses 
by  Hooke's  Law. 

^x  -  ^S  =  ^^"'X  ^  \   ^^X^ 
(2.5)         Oy  -  v,f^  =  E[V,Y  +  \   W?y] 

^xy  =  2(fTvT  ^^^Y  +  ^'X  "^  ^'X^'Y^  ' 
The  bending  moments  are  determined  from  W(X,Y)  by  the  formulas, 

M^  =  -tT[W,j^  +  vW,^] 
(2.^1)  My  =  -tT[W,,^Y  +  vW,^^] 

\y=  (l-v)tTW,xY  • 

It  is  assumed  that  the  boundary  of  the  plate  Is  simply 
supported.  Then  (2.4)  imply  that, 

(2.5a)  W  =  AW  =  0  on  S  . 
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It  is  further  assumed  that, 

^  vTY 

cr^  =  -  T,   Y  =  -~,    for  X  =  0,a,  0  l  Y  <  b  ; 

(2.5b) 

TX 
<Jy  =  0,  "■]]■=  -    g-,  for  0  1  X  <  a,  Y  =  0,b  . 

The  bounlary  conditions  are  especially  convenient,  although 
not  essential,  for  the  analysis  given  in  the  next  section  anl 
the  numerical  computations  lescribe:!  in  the  last  three  sections 
of  this  paper. 

A  solution  of  the  boundary  value  problem  (2.1)  (2.5) 
valid  i:or  all  values  of  the  parameters  a,  b  an1  T  is  given  by 
F°(X,Y),W°{X,Y),  U°(X,Y)  and  V°{X,Y)  where, 

i     (2.6)    p°  ^.  .  'L^l,    w°  ^  0,  U°  ^  -  TX^  v°  -  ^^ 

This  is  called  the  "unbuckle i"  solution  since  the  mi  Isurface 
of  the  plate  remains  in  the  X,Y  plane.   To  determine  other 
solutions,  it  is  convenient  to  define  the  following  dimension., 
less  variables  and  parameters: 

X  =  Xb"^,  y  =  Yb-\  £   =   a/b,  c  =  [12( l-v^) ]^/2 

u(x,y)  H  c2bt-2[Tj(x,Y)-U°{X)];  Vtjc.y)  =  c%t-^[V{X,Y)-V°  {Y)] , 
(2.7)  ,  2 

W(x,y)  =  ct"  M(X,Y),  r(x,y)  =^  [P(X,Y) -P°(X,Y) ] , 
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Thus  the  iimensionlcss  indepenclent  variables  are  in  the  ranges, 

0  ±  X  ±  £ ,    0:ly_il  where  &   Is  the  aspect  ratio  ox"  the 

plate.   The  jarameter  A  is  called  the  load.   In  terms  of  the 

"excess"  stress  function  r(x,y)   the  membrane  stresses  and  the 

corresponding  reduced  stresses  )       ,  )        and  )     are  given  by 

' — X  — y     ' — xy 

2      2 

YI^(x,y)  =  ^(|)(r^(X,Y)  =   -   -K   +   f.  ^, 
— ^A       g  t  ■'^  yj 

2   2 

(2-8)  Hyt'^.y)  ^  |-(|)   o-yix.Y)  =  f,_^^, 

Using  (2,7)  and  (2.8),  Hooke's  lav;  (2.3)  reduces  to 

1   2 
f,    -  vf,    =  u,   +  —  w.    , 
'yy     'xx    'x   2   ^ 


2 


(2.9)         f,    -  vf,    =  V,   +  -  w, 

'xx      'yy     y        2        y 

f,    =  -  [S( 1+v) ]""[u,  +v,  +w,  w,  ]  . 
'  xy     u  V  J.  V  /  J   '■  '  y   ^  X   X  '  y  ■■ 

The  reduce:!  moments  are  given  in  terms  of  the  reduced  lateral 
displacement  by 
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3,  2 
my(x,y)  ^  ^^   My(X,Y)  =  w,^^  +  vw,^^  , 

Et 

3  2 
m^^(x,y)  =  -^-5 ^  M   (X,Y)  =  w      . 

In  terms  of  the  dlmensionless  variables  (2.7)  the  von 
Karman  equations  (2.1)  become 

(2.11a)  ^  A  V/+AW,   =  x"",  w,  +f,  w,  -2r,  w,   =  H[w,r]  , 

^      '         'xx    'yy  'xx   'xx  'yy    'xy  'xy     l  ,  j  , 
(2.11b)    A^r  =  w?   -  w,   w,    =   K[vj]  , 

where  A  Is  the  Laolacian  with  respect  to  x  and  y.   Employing 
(2.7)  an-1  (2.3),  the  boundary  conditions  (2.5)  become 

w  =  Aw  =  0  on  S, 

-A   ,  V  r=  0,   for  x  =  0,  £,    O^y^l  , 


(2.12) 

■ 

^     :. 

(2.13) 


>    =  0,  u  =  0,   for  01x<^,  y=0,  1  . 

y 


Using  (2.8)  and  (2.9)  the  conditions  (2.12)  and  (2.13)  imply 


(2.14)         f,^  =  ^^yy  =  °  °^  ^  . 
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The  conditions  (2.l4)  can  be  further  simplified  to  yield 

(2.15)  f  =  Af  =  0  on  S   . 

To  show  this,  the  tangential  derivatives  in  (2.l4)  are  inte- 
grated along  each  edge  of  S.   Pour  of  the  eight  constants  of 
integration  are  eliminated  by  assuming  that  f(x,y)  is  continuous 
on  S.   Since  f(x,y)  is  determined  from  (2.11-2.13)  to  ^^;ithin 
an  arbitrary  linear  function  of  x  and  y,  three  of  the  remaining 
constants  are  eliminated.  Thus  v\fe  obtain,  for  example,  that 
on  S, 

(2.16)  f(x,0)  =  f(^,y)  =  0,  f(0,y)  =  Ay,  f(x,l)  =  -  |(x-f,), 

where  A  is  an  arbitrary  constant.   Since  it  is  required  that 
every  solution  of  the  boundary  value  problem,  (2.11),  (2.12), 
(2.14)  and  (2.l6)  be  a  solution  of  the  boundary  value  problem 

(2.11)  -  (2.13),   A  must  be  zero.   Otherwise,  w  =  0  , 
f(x,y)  =  -  ^(x-J3)y  is  a  solution  of  (2.11),  (2.12),  (2.l4) 
and  (2.l6)  that  violates  the  conditions  on  u  and  v  in  (2.13). 
Thus  the  conditions  (2=15)  follow  from  (2.l4)  and  (2.l6).   The 
differential  equations  (2.11)  and  the  boundary  conditions 

(2.12)  and  (2.15)  are  used  in  the  remainder  of  the  paper.   The 
displacements  u  and  v  are  determined  from  a  solution  of  this 
problem  by  Integrating  (2.9)  and  using  the  conditions  on  u  and 
V  in  (2.13). 
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It  can  be  shown  using  (2.12)  and  (2.15),  that  the 
difference  between  the  potential  energies  of  a  buckled  and 
the  unbuckled  state  is  proportional  to  the  functional  Q, 

where 

1  H 

(2.17)      Q  ^  y  y^  [(Aw)^-Aw?^  +  (Af)^]dxdy  . 

o  o 
Q  is  hereafter  referred  to  as  the  energy. 

3.   The  Unbuckled  State. 


In  this  section  we  prove  that  if  A  is  less  than  or 
equal  to  the  minimum  eigenvalue,  'hi^) ,    of  the  classical 
linearized  buchling  theory  [l],  then  the  nonlinear  boundary 
value  problem  has  the  unique  solution  f  =  w  =  0  .   The 
linearized  theory  is  obtained  by  setting  H  =  K^  0  in  (2.11). 
Then  f  =  0  ,  and  w  is  a  solution  of  the  eigenvalue  problem: 

2 

(5.1)  A  w  +  w,    =0,   v/  =  Aw  =  0  on  S  . 

This  problem  has  non- trivial  solutions  ( eigenf unctions ) , 
(5.0)   w  =  Wj^(^'y)  ^  Aj^^sin(mT7x/^)  sin  njry,   m,n  =  1,2,... 
if  and  only  if  ■:  ■   . 

(3.2)  ^  =  ^mn  -  (V-e)^[m  +  n^^^/m]^  ,  m,n  =  1,2,  ... 


Here  A   are  arbitrary  constants  and  A^^^  are  the  eigenvalues, 
Thus  it  follows  from  (3.2)  that 
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(3.3)  Mi)    =   mm.  -K^^      . 

m 

We  can  prove,  employing  classical  methods  [12]  and  (3.1)- 

(3.3),  that  the  inequality 

1  I 
(3.^)         f  !    [(Aw)^  ->,w?^Jdxdy  >  0 

o   o 
holds  for  all  "admissible" functions  w  ^  0  if  A  <  A*   ^^  also 

holds  for  X  =  a  provided  that  w  ^  0,  v/   .   If  w  =  w   and 
—  ~   '   mn       ~  mn 

A  =  _X,  (3.^)  reduces  to  an  equality. 

To  prove  the  uniqueness  theorem,  the  boundary  value  problem 

B,  is  formulated  as  follows:  Find  an"admissible"  function  w(x,y) 

which  satisfies  the  differential  equation  (2.11a)  and  the 

boundary  conditions  (2.12).   In  (2.11a)  and  (2.17)  f(x,y)  is 

considered  as  a  functional  of  w(x,y)  defined  by  (2.11b)  and 

(2.15).   Hence  Q  is  a  functional  of  vj(x,y)  only.   It  can 

easily  be  shovm  that  if  v;  is  a  solution  of  B,  Q  is  stationary, 

i.e.  the  first  variation, 6Q,  of  Q  vanishes.   Furthermore,  if 

w  is  a  solution  of  B  it  is  also  an  admissible  variation.   Hence, 

"   setting  the  variation  of  w  equal  to  w,  the  equation 

5Q  =  0  becomes, 
1  £ 
(3.5)      r  /[(^")^  -  ^w?.^  +  2(Af)^]dxdy  =  0  . 
-       00 
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Thus  it  follows  from  {^A)   and  (3-5)  that  if  A  <  A>  w(x,y)  =  0 

and  hence  f (x,y)  =  0  .   If  A  =  A^  then  either  w(x,y)  =  0  or  it 

is  an  eigenfunction  of  (5.1).   Suppose  that  A  =  A  and  w(x,y)  = 

w   (x,y) .   Then,  from  (5.5)  >  the  remarks  follov/lng  (5.'^)  and 

from  (2.11b)  and  (2.15)  it  follows  that  w   satisfies 

mn 

K[v;   ]  =  0  for  all  x  and  y  in  D,   This  is  impossible  and  hence 
mn  "^ 

the  proof  is  complete ^ 

The  above  proof  is  similar  to  one  given  by  Friedrichs 
and  Stoker  [2]  in  their  study  of  the  symmetric  bucicling  of 
circular  plates. 

4.   Numerical  Methods. 

To  obtain  approximate  solutions  oC  the  boundary  value 
problem  it  is  convenient  to  first  transform  the  two  fourth 
order  differential  equations  (2.11)  into  a  system  of  four 
second  order  equations.   Thus  the  functions  0(x,y)  and 
(j)(x,y)  are  defined  by 

(4.1)  C(x,y)  =  Aw(x,y)  ,  (|)(x,y)  =  Af(x,y)  . 

Hence  (2.11)  is  written  as,  :.:■.  '■ 

A(|)  =  K[w]  ,  Af  =  (|) 


(4.2) 


and 


An  =  -Aw,  ^  +  H[w,f],     Aw  =  O  , 
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(4.3)  (l)  =  f  =  O=w  =  0onS. 

Starting  from  an  initial  estimate  v/^°'(x,y)  of  the 
solution  for  fixed  A  and  £      ,    approximate  solutions 
of  (4.2)  and  (4.3)  are  sought  by  an  "accelerated"  iteration 
procedure.   A  sequence  of  Iterates  [w^'^'(x,y),  f^'^'(x,y)  , 
^        (x,y),0^   (x,y)]  is  defined  by  the  recursions, 

A(t)^^)  =  KLw^'^h  ,  (|)^^^  =  0  on  S  , 

Af(")  =  c|)("^  ,  f(^)  =  0  on  S  , 

(4.4)  AO^"^  =  -Aw!;i  +  H[w(^^f(^h,   o(")  =  0  on  S  , 

_(n+l)    (^  _(n) 

Aw     =  O^  '  ,  w    =  0  on  S  , 

where  6  is  the  acceleration  parameter.   If  the  Initial  iterate 
w   (x,y;,  is  sufficiently  simple,  say  a  polynomial  in  x  and 
y,  it  is  conceivable  that  each  successive  Itei'ate  can  be  ex- 
plicitly obtained.   However,  in  general,  they  will  rapidly 
become  unwieldy  and  unsuitable  for  numerical  tabulation. 
Therefore,  the  iterates  are  computed  approximately  by  a  finite 

difference  procedure. 

Each  iterate  in  (4.4)  is  determined  as  the  solution  of 
a  linear  boundary  value  problem  of  the  form, 

(4.5)   ■  ,,   Ag(x,y)  =  c(x,y)  ,    g  =  0  on  S   .     . 

Here  c(x,y)  is  a  known  function.   To  apply  the  finite  difference 
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procedure  to  (4.5),  the  region  D  with  boundary  S  Is  covered  by 
a  uniform  rectilinear  mesh  with  spacing  6  in  the  x  and  y 
directions.   The  mesh  is  selects i  so  that  the  boundary  lines 

of  D  coincide  with  mesh  lines.   There  are  M  +  1  and  N  +  1 

lines  in  the  x  and  y  directions  so  that 

(4.6)  5  =  VM  =  1/N 

The  mesh  points  (x.,y.)  are  defined  as, 

(4.7)  x^  =  15  ,  y.  =  j5  ,  i  =  0,1, ...,M,  j  =  0,1, ...,N  . 

The  set  of  mesh  points  which  are  interior  points  of  D  is 
called  the  mesh  region  D^  .   The  set  of  the  mesh  points  on 
S  is  denoted  by  S^^.   A  function  defined  only  at  the  points  of 
Dc-  is  called  a  mesh  function.   It  is  assumed  that  at  each  point 

0 

of  Do  the  solution  of  (4.5),  g(x.,y.),  is  approximated  by  the 
mesh  function  'IS-.t't  which  satisfies  the  difference  equations 

Vij  =  ^ij  '  (^i.yj)^D6    ' 

(4.8) 

g^j  =  0  (x,,yj)tS5   . 

Here  A^-  is  the  nine  point  difference  approximation*  [13]  to 
the  Laplacian  and  is  defined  by 


*  Numerical  experiments  conducted  by  S.  Schechter  (privately 
communicated  to  the  authors)  for  the  solution  of  Laplace's 
equation  indicate  that  the  nine  point  scheme  yields,  roughly, 
the  same  accuracy  as  the  conventional  five  point  scheme  with 
twice  as  many  mesh  points. 
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c 
The  mesh  functions  ^;  c .  .  \    are  obtained  in  each  case  by 

replacing  the  second  derivatives   occurring  in  the  right 

hand  sides  of  (4.4)  by  their  centered  second  difference 

approximations . 

The  difference  equations  (4.S)  and  (4.9)  yield  a  system 

of  (M-1)  X  (N-1)  linear  algebraic  equations  for  the  mesh 

function  <g.  .\        •      The  solutions  of  this  system  converge 

to  the  solution  of  (4.5)  as  a  B  -^  0  .   Hence,  for  a  sufficiently 

small  5  we  expect  that  the  solution  of  the  algebraic  problem 

yields  a  "close"  approximation  to  the  solution  of  the  boundary 

value  problem  (4.5).   The  coefficient  matrix  of  the  algebraic 

system  is  triaiagonal  with  respect  to  matrices,  or  of  quasi- 

tridiagonal  form.   This  system  is  solved  by  a  direct  method 

[l4]  in  which  the  coefficient  matrix  is  factored  into  the 

product  of  two  matrices,  one  lower  and  one  upper  triangular 

with  respect  to  matrices.   The  solution  is  obtained  by 

successively  inverting  the  triangular  systems.   This  requires 

the  inverses  of  N-1,  (M-l)  by  (M-1)  matrices.   They  are  obtained 

by  Gauss  elimination  vjith  pivotal  condensation.   For  a  fixed 

6,  the  iterative  procedure  (4.4)  requires  that  these  inverses 

be  computed  only  once.   Thus  for  each  iteration  except  the 
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rirst,  only  multiplication  of  vectors  by  matrices  and  evalua- 
tions of  th2  inhomogeneous  terms  in  (4.4)  are  needed  to 
compute  an  Itei^ate. 

Wl:ier.  the  iterations  have  converged  for  specified  values 
of  "h     and  £,    the  stresses  and  the  potential  energy  are 
computed  from  finite  difference  equivalents  of  (2.8),  (2.10) 
and  (2,].7).   The  integrals  in  (2.17)  are  approximated  by 
Simpson's  rule. 

5^   Computational  Procedures. 

All  numerical  computations  were  performed  on  the  IBM  7090 
computer  at  the  Courant  Institute   of  Mathematical  Sciences. 
The  N-1  matrices  required  for  the  solution  of  (4.8)  and  (4.9) 
are  stored  in  the  internal  or  fast  access  memory.   One 
iteration,  i.e.  a  computation  of  w    '    from  w^   ,  requires 
approximately  15  seconds  for  a  27  by  27  mesh.   The  25  inver- 
sions require  approximately  1.2  minutes.  The  internal  storage 
of  the  N-l  inverted  matrices  limits  the  size  of  the  mesh. 
Thus  for  the  square  plate  we  are  limited  to  625  interior  mesh 
points  and  for  a  rectangular  plate  with  &  =  2     to  74l  points. 
Finer  meshes  can  be  employed  if  tape  storage  is  used.   However, 
this  considerably  increases  the  computing  time. 

Application  of  the  difference  equations  (4>,8)  and  (4.9) 
requires  the  mesh  spacing  to  be  "sufficiently  small". 
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A  satisfactory  value  of  5   is  determined  by  comparing  the 

solutions  of  a  sequence  oT  test  calculations  witli  successively 

finer  meshes.   VJe  conclude  from  the  results  of  these  tests 

that  the  625  point  mesh  (5  =  — ),  and  the  7^1-1  point  mesh 

26 

(6  =  1/20)  provide  sufficiently  accurate  solutions  for  plates 
with  aspect  ratios  of  £  =  1  and  i  =  2. 

The  following  numerical  convergence  criterion  is  employed: 

/c  1  \  I  (n+1)   (n) I  . 

(5-1)  max     |w^  .   '-w:  .' |  <  e  , 

(x.,yj)tD  ^J     ^^ 

where  e  =^  0  is  a  prescribed  "small"  number.   For  all  the 

_7 
calculations  epsilons  were  in  the  range  1.0xlO:le:Sl.7x 

_5 
10   .   In  general,  permissible  values  of  e  increase  with  A. 

Starting  from  a  given  initial  iterate,  the  number  of 

iterations  needed  to  satisfy  (5«1)   depends  on  the  value  of  Q 

used  in  (4.4).   The  optimal  value,  6    ,    is  defined  as  that 

value  which  minimizes  the  number  of  iterations.   To  obtain 

estimates  for  6  we  resort  to  a  series  of  test  calculations. 

o 

In  general,  we  find  that  the  estimated  Q  decreases  as  A 
increases  and  it  also  depends  on  which  of  the  non-unique 
solutions  is  being  approximated. 

The  initial  iterate  should  be  as  close  an  approximation 
to  the  solution  as  possible.   To  obtain  w    for  the  bifurcating 
solutions  and  for  a  fixed  value  of  i  and  any  A  we  proceed  as 
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follows:   For  A  :!  A,w  =  0  Is  the  only  solution,  thus  only 

A  >  A  is  considered.   A  value  of  A  is  chosen  which  is 

"slightly  greater"  than  A   and  w^°'  is  assumed  in  the  form 

of  a  corresponding  eigenfunction.   Then  the  iterations  are 

performed  until  they  converge  to  a  buckled  solution.   This 

solution  is  then  used  as  the  inital  iterate  for  the  load 

A  +  AA  where  AA  is  a  "small"  increment.   In  this  manner 

solutions  are  obtained  for  increasing  and  decreasing  sequences 

of  loads  on  a  specific  branch.   As  A  -^  A    the  number  of 

mn 

iterations  necessary  for  convergence  rapidly  increases.   In 

fact,  for  A  sufficiently  "close"  to  A   we  are  unable  to 

''  mn 

obtain  convergence.* 

For  large  loads  the  number  of  iterations  increase  as 

A  increases.   We  also  find  that  there  are  A's  at  which  the 

iterations  jump  or  converge  to  a  solution  which  bifurcates 

from  another  eigenvalue.   This  indicates  the  original  solution 

no  longer  exists  at  that  value  of  A,  vjhich  is  unlikely,  or  that 

the  second  solution  is  more  stable  with  respect  to  our  iterative 

procedure.   The  results  show  that  the  jumping  occurs  only  when 

the  energy  of  the  second  solution  is  less  than  that  of  the 

first.   In  one  case,  discussed  in  the  next  section,  the 

iterations  jump  to  an  asymmetric  solution. 

*In  this  range  approximations  can  be  determined  for  solutions 
bifurcating  from  simple  eigenvalues  by  a  perturbation  procedure. 
The  results  of  this  calculation  are  not  discussed  in  this 
paper. 
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6.   Presentation  and  Discussion  of  Results*. 

The  numerical  method  discussed  in  the  two  preceding 
sections  is  applied  to  the  square  plate,  i  =  1,  and  to  the 
rectangular  plate  with  i  =  2.      For  the   square  plate  we 

consider  the  solutions  which  apparently  bifurcate  from  the  tv;o 

2  ^S  2 
lowest  eigenvalues,  A  =  47r  ,  -7-f^7r   .   They  are  called,  respective- 
ly, the  first  and  second  modes.   Numerical  solutions  are  ob- 
tained for  loads  which  are  slightly  greater  than  the  eigen- 
values and  ±   700  .   For  39.5  ^  A  <  76  the  Iterations  always 
converge  to  the  first  mode. 

For  the  rectangular  plate  approximations  are  obtained 
to  the  asymmetric  solution  and  to  the  solutions  which  seem 
to  bifurcate  from  the  four  lowest  eigenvalues 

A  =  ^TT^,  -^TT^,  ^TT^,   ---TT^  .  'K   =  "Ptt^  Is  a  double  eigenvalue 
and  is  obtained  from  (3.2)  with  m  =  1,  n  =  1  and  m  =  4,  n  =  1. 
The  solutions  corresponding  to  the  first,  second  and  fourth 
eigenvalues  are  denoted  as  the  first,  second  and  fifth  modes 
and  the  solutions  corresponding  to  the  third  eigenvalue  with 
m  =  4  and  m  =  1,  as  the  third  and  fourth  modes  respectively. 

The  range  of  loads  for  which  solutions  have  been  com- 
puted for  the  rectangular  plate  are  shown  in  Table  I.   Using 
the  second  mode  solution  for  A  =  52  as  the  initial  iterate, 
the  Iterations  converge  to  the  first  mode  at  A  =  51  •   Similar 
transitions  are  observed  for  other  modes. 

*  All  calculations  employed  v  =.3.   The  authors  acknowledge 
the  assistance  of  C.  Kaman  and  M.  Greenberg  in  reducing  the 
data. 
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They  are  indicated  in  the  remarks  column  of  Table  I. 

Table  I 
Limits  of  A  for  computed  solutions  for  the  rectangular  plate, 


Mode 

Lower  Limit 

1 
Upper  Limit 

Remarks 

of  A 

of  A 

1 

1 

59.5       \ 

120 

at  A  -. 
at  A  = 

=  39-0  -^  zero 
=  120  -?  A-Mode 

A 

110 

400 

at  A  = 

.  450  -^  3rd  Mode 

2 

52 

500 

at  A  = 

=  51   ->  1st  Mode 

5 

112 

500 

4 

125 

390 

at  A  = 
at  A  = 

=  120  -^  2nd  Mode 
-   400  -^  5th  Mode 

5 

550 

500 

at  A  = 

=  300  -^  2nd  Mode 

In  each  case  the  transition  occurs  to  a  state  of  equal  or  lov/er 
energy.   For  A  in  the  range  110  ^  A  i  120  both  the  first  and  the 
asym.metric  modes  are  determined.   Although  their  energies  are 
essentially  equal  , the  shapes  differ ' somewhat.   The  asymmetric 
modes  appears  to  bifurcate  from  the  first  mode  below 
A  =  110  . 

In  Pig.  1  the  variation  v/ith  A  of  the  energies  for  the 
various  modes  is  shown.   Since  the  values  of  these  energies 
are  relatively  close  to  each  other,  we  prefer  to  shov;  in  Fig's, 
lb  and  Ic  the  energies  relative  to  the  second  mode. 
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The  second  mode  energies  for  the  square  and  rectangular  plates 
are  illustrated  in  Pig.  la. 

At  A  ^  130,  550  the  first  and  second  modes  of  the  square 
plate  have  equal  energies.   Similarly  for  the  rectangular 
plate,  the  first  and  second  modes  have  equal  energies  at 
7\   ^■:  96,    the  asymmetric  and  third  modes  have  equal  energies  at 
7\  ^  248.  We  now  employ  these  observations  and  describe  a 
possible  mechanism  to  explain  the  experimentally  observed  mode 
jumping  [7]. 

The  same  mechanism  was  proposed  in  [8]  to  explain  the 
snap  buckling  of  spherical  caps.   It  is  essentially  a  generali- 
zation of  Priedrlchs  [15]  concept  of  an  intermediate  buckling 
load  for  spherical  shells.   A  buckled  state  B. (A)  with  energy 
Q.  (A)  is  said  to  be  preferred  to  a  state  B,  (?v)  with  energy 
Q^(A)  if  0^(7\)  <  Q.^(A).   The  "energy-crossing"  loads  are 
defined  as  those  loads  for  which  Q.  =  Q,  .   V/e  assume  that  it 
is  always  possible  for  the  plate  to  jump  from  any  state  to  a 
preferred  one.   If  state  B. (A)  is  preferred  to  B,  (A)  and  they 
are  "close",  i.e.  they  have   similar  stresses  and  displacements, 
then  "small"  disturbances,  which  are  always  present  in  an 
experiment,  could  initiate  the  Jumping  from  B,  (A)  to  B. (A)  . 
If  B.(A)  and  B,  (A)  are  not  close  and  a  sufficiently  large 
disturbance  is  not  present,  it  may  still  be  possible  to  jump 
from  B,  (A)  to  B.(A)  via  an  intermediate  state  B.(A). 
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Thus  If  Q.(A)  <  Q,(A)  <  0.   (A)   and  B.(A)   is  close  to  B,  (A) 
then  a  small  disturbance  could  initiate  a  jump  from  B,  (A)  to 
B.(A).   Once  the  plate  is  in  motion  it  may  then  be  possible 
to  reach  the  non-adjacent  preferred  state  B.(A).   The  asymet- 
ric  solution,  which  is  found  for  the  rectangular  plate,  could 
conceivably  serve  as  an  intermediate  state  between  the  first 
mode  and,  say,  the  third  mode.   To  Illustrate  this,   w(x,l/2) 
is  shovm  in  Pig.  2  for  the  first  mode  at  A  =  120,  for  the  asym- 
metric mode  for  a  sequence  of  loads,  and  for  the  third  mode  at 
A  =  400. 

If  this  is  the  true  jumpting  mechanism  then  the  energy- 
crossing  loads  provide  lower  bounds  for  the  experimentally 
observed  jumping  loads.   If  tne  states  B.(A)  and  B.  (A)  are 
not  close  or  if  there  are  no  close  intermediate  states  then 
jumping  becomes  more  difficult  and  the  crossing  loads  may  give 
inaccurate  bounds. 

Some  of  the  results  obtained  for  the  square  plate  with  A 
in  the  range  4o  ^  A  j^  700  for  the  first  mode  and  in  the  range 
8o  jS  A  :!  700  for  the  second  mode  are  now  briefly  described. 
For  A  <   59.5  the  iterations  always  converged  to  the  unbuckled 
solution. 

Representative  graphs  of  the  variation  with  x  at 
y  =  1/2  of  the  dimensionless  stresses  and  displacements  are 
shown  in  Pig='s  3-^6. 
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The  first  (second)  mode  Is  symmetric  ( ant i- symmetric)  with  respect 

to  the  line  x  =  1/2  .   Hence,  only  one-half  of  the  plate 

0  ^  X  jf^  1/2  is  shown.   The  formation  of  boundary  layers  as 

>v  Increases  is  clearly  evident  in  these  figures.   The  maximum 

of  the  absolute  values  of  each  of  the  quantities  shown  moves 

towards  the  edge  as  A  increases. 

The  dimensionless  bending  stresses  m   (not  shovm)  are 

similar  in  shape  to  m  but  considerable  smaller  in  magnitude. 

The  stresses,  }        and  Z       ,  of  the  second  mode  are  similar  to 
'  X       ^y 

those  of  the  first  mode,  Pig's.  5  and  6  and  are  therefore  not 
shown.   To  further  illustrate  the  foimation  of  boundary  layers 
near  the  loaded  edges  we  consider  the  energy,  Q(x;7\),  along 
lines  X  =  constant  which  is  defined  by  the  integral* 


J 


Q(x;A)  =    r[(Aw)^  -  Tvw,^  +  (Af)^]dy 


Approximations,  Q,(x,  ;A)  to  Q(xjA)  are  obtained  along  each  mesh 
line  parallel  to  the  loaded  edge  by  approximating  the  integral 
using  Simpsons  rule  and  evaluating  the  derivatives  with 
centered  difference  quotients.   The  variation  of  Q(x.;a) 
with  X  is  sho'iA?n  in  Fig.  7  for  the  first  mode  for  a  sequence 
of  X  values.   Corresponding  graphs  (not  shown)  for  the  second 


*The  idea  of  examining  the  variation  of  Q(x;A)  was  essentially 
suggested  to  the  authors  by  J.J.  Stolcer. 
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mode  of  the  square  plate  and  the  six  solutions  computed  for 
the  rectangular  plate  are  similar  to  those  shown  in  Fig.  7. 
The  numerical  results  indicate  that  for  all  solutions  considered 
the  variation  of 


max    Q(x. ;A)  =  Q(0;A) 

0  <  1  <  N-1    ^ 

with  A  can,  for  sufficiently  large  'K,   be  approxim.ated  by, 

Q(0;A)  ^  }J^^/Q0  /2      . 

This  formula  may  provide  a  useful  basis  for  estimates  of  the 

maximum  load  that  the  plate  can  support j  see  also  the  discussion 

below  of  von  Karmans  ultimate  load  conjecture. 

The  variation  with  y  of  the  energy  along  lines  parallel 

to  the  loading  direction  gives  no  indication  of  boundary  layers 

near  y  =  0,  and  y  =  1  for  the  range  of  loads  that  are  considered, 

In  fact  the  boundary  layers  near  the  unloaded  edges  are  much 

"weaker"  and  less  developed  than  the  ones  near  the  loaded  edges. 

For  example,  there  appears  to  be  no  boundary  layer  development 

for  T^  near  y  =  0  and  only  a  weak  one  for  the  first  mode  >~ 

as  shown  in  Fig.  8  and  none  for  the  second  mode  ~  (not  shown) . 

However,  for  the  dimensionless  shear  stress  y~       an  "internal" 

-^^xy 

boundary  layer  seems  to  develop  at  x  =  1/2,  see  Fig.  9« 

The  variety  of  solutions  and  the  mode  jumping  indicate 
that  the  asymptotic  description  as  A  ->  oo  of  the  solutions  of 
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the  nonlinear  boundary  value  problem  may  be  quite  complicated. 
Von  Karman  [l6]  proposed  an  analysis  for  the  buckled  plate  using 
boundary  conditions  other  than  (2.15).   In  this  analysis  the 
plate  is  deformed  by  a  uniform  normal  membrane  displacement. 
It  is  assumed  that  for  displacements  considerably  in  excess  of 
the  critical  one,  the  load  is  carried  by  two  narrow  strips 
parallel  and  adjacent  to  the  unloaded  edges.   The  load  is  also 
assumed  to  be  uniformly  distributed  across  the  strips.   The 
graphs  in  Fig.  10  indicate  that  von  Karman' s  conjecture  miay 
be  reasonable  away  from  the  loaded  edges  for  the  boundary  con- 
ditions considered  in  this  paper.   The  corresponding  graphs 
for  the  second  mode  and  for  the  rectangular  plate  have, 
roughly,  the  same  form  as  in  Pig.  10.   Thus  the  von  Karman 
formula  [l6]  may  also  provide  an  approximation  for  the  ultimate 
load  for  the  plate  problem  investigated  in  this  paper. 

The  second,  fourth  and  fifth  (first  and  third)  modes 
for  the  rectangular  plate  £  =  2  are  symmetrical  (antisymmetrlcal) 
with  respect  to  the  line  x  =  1/2  .   For  half  the  plate, 
0  ±  X  <   1/2  ,  and  the  range  of  loadings  considered  the  stresses 
are  similar  in  form  to  each  other  and  to  those  of  the  square 
plate.   In  Fig's.  11-13,  representative  results  are  shown  for 
the  rectangular  plate  with  A  =.  350. 
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mode  of  the  square  plate  for  an  increasing  sequence 

of  A. 
Pig.  4.    The  variation  of  m  (x,l/2)  with  x  for  the  first 

mode  of  the  square  plate  for  an  increasing  sequence 

of  A. 
Pig.  5.    The  variation  of  7~  (x,l/2)  with  x  for  the  first 

mode  of  the  square  plate  for  an  Increasing  sequence 

of  A. 
Pig.  6.    The  variation  of  7~^(x,l/2)  with  x  for  the  first 

mode  of  the  square  plate  for  an  increasing  sequence 

of  A. 
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Pig.  7-   The  percentage  of  the  total  energy  that  is  contained 

in  strips  perpendicular  to  the  loading  direction. 
Pig.  8.   The  variation  of  >~  (l/2,y)  vjith  y  for  the  first 

mode  of  the  square  plate  for  an  increasing  sequence 
of  A. 
Pig.  9.   The  variation  of  ^^I   (7/26, y)  with  y  for  the  first 
mode  of  the  square  plate  for  an  increasing  sequence 
of  -\. 
Pig.  10.  The  variation  of  7"  {x,y)  with  y  for  x  =  1/26,  7/26, 
1/2  for  the  first  mode  of  the  square  plate  for 
A  =  500. 
Pig.  11.   The  variation  of  w(x,l/2)  with  x  for  the  rectangular 

plate,  &   -T  2,  for  A  =  350. 
Pig.  12.   The  variation  of  m  (x,l/2)  with  x  for  the  rectangular 

plate,  A   =  2,    for  A  =  350. 
Pig.  13.   The  variation  of  y~  (x,l/2)  with  x  for  the  rectangular 
plate,  £  =  2,  for  A  =  350. 
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